Introduction
Tutorial for calculating and plotting dNBR Fire Severity Index using Landsat 8/9 satellite imagery.
Additional Contributors
Conservation Data Lab, Randy Swaty, Sam Ettenborough
Setup
Please review the README.md document for installation and setup instructions.
Importing Landsat Files
Code
# GEE project title
project <- "ee-samettenborough"
# fire 1 coordinates
fire_1 <- list (- 91.253 , 16.8194 , - 89.4018 , 18.2112 )
# these two fires started end of may and last measurements were in end of july
aoi_03 <- here:: here ("data/areas-of-interest/EMSR727_AOI03_BLP_PRODUCT_areaOfInterestA_v1.shp" )
aoi_05 <- here:: here ("data/areas-of-interest/EMSR727_AOI05_BLP_PRODUCT_areaOfInterestA_v1.shp" )
# pre fire start and end date search window
pre_start_date <- "2024-01-01"
pre_end_date <- "2024-03-01"
# output folder
out_dir <- "data/pre"
# uncomment this section to re-download the images
# process_aoi(
# project,
# pre_start_date,
# pre_end_date,
# out_dir,
# geometry = aoi_05,
# prefix="aoi_05_landsat_pre_SR5_SR7_",
# )
# pre fire start and end date search window
post_start_date <- "2024-04-01"
post_end_date <- "2024-08-01"
# output folder
out_dir <- "data/post"
# # uncomment this section to re-download the images
# process_aoi(
# project,
# post_start_date,
# post_end_date,
# out_dir,
# geometry = aoi_05,
# filter_lim = 30,
# prefix="aoi_05_landsat_post_SR5_SR7_",
# )
Visualize Pre and Post Images
Code
pre_fire_files <- list.files ("data/pre/" ,
pattern = "aoi_05_landsat_pre_SR.*" ,
full.names = TRUE )
pre_fires <- pre_fire_files |>
map (rast) |>
reduce (mosaic)
# Plotting in RGB to test
par (col.axis = 'white' , col.lab = 'white' , tck = 0 )
plotRGB (pre_fires,
r= 2 , g= 1 , b= 1 ,
stretch = 'lin' ,
axes = TRUE ,
main = 'Prefire RGB Composite Image, Bands 5,7' )
Code
# # Create leaflet map
# leaflet() %>%
# addTiles() %>%
# addRasterImage(
# pre_fires,
# opacity = 0.7,
# group = "Pre-Fire RGB Composite"
# )
Code
post_fire_files <- list.files ("data/post" ,
pattern = "aoi_05_landsat_post_SR.*" ,
full.names = TRUE )
post_fires <- post_fire_files |>
map (rast) |>
reduce (mosaic)
par (col.axis = 'white' , col.lab = 'white' , tck = 0 )
plotRGB (post_fires,
r= 2 , g= 1 , b= 1 ,
stretch = 'lin' ,
axes = TRUE ,
main = 'Postfire RGB Composite Image, Bands 5,7' )
Calculate Normalized Burn Ratio and Burn Severity (dNBR)
Code
# Calculating the Delta Normalized Burn Ratio Index
nbr_pre <- 1000 * (pre_fires[[1 ]] - pre_fires[[2 ]]) /
(pre_fires[[1 ]] + pre_fires[[2 ]])
nbr_post <- 1000 * (post_fires[[1 ]] - post_fires[[2 ]]) /
(post_fires[[1 ]] + post_fires[[2 ]])
# Attempt to calculate the Differenced Normalized Burn Ratio Index
dnbr <- (nbr_pre)- (nbr_post)
Display Pre and Post Fires
Code
nbr_stack <- c (nbr_pre, nbr_post)
names (nbr_stack) <- c ("Pre-fire NBR" , "Post-fire NBR" )
nbr_stack_df <- rasterdf (nbr_stack)
ggplot (nbr_stack_df) +
geom_raster (aes (x = x,
y = y,
fill = value)) +
scale_fill_gradient (name = "NBR" ,
low = "lightyellow" ,
high = "darkgreen" ) +
coord_sf (expand = FALSE ) +
annotation_scale (location = 'bl' ) +
facet_wrap (facets = vars (variable),
ncol = 1 ) +
theme_void ()
Display dNBR
Code
dnbr_df <- rasterdf (dnbr)
ggplot (dnbr_df) +
geom_raster (aes (x = x,
y = y,
fill = value)) +
scale_fill_gradient2 (name = "DNBR" ,
low = "blue" ,
high = "red" ,
midpoint = 0 ) +
coord_sf (expand = F) +
annotation_scale (location = 'bl' ) +
theme_void ()
Classify Severity
Code
# Classifying the index values into a matrix
rclas <- matrix (c (- Inf , - 970 , NA , # Missing data
- 970 , - 100 , 5 , # Increased greenness
- 100 , 80 , 1 , # Unburned
80 , 265 , 2 , # Low severity
265 , 490 , 3 , # Moderate severity
490 , Inf , 4 ), # High severity
ncol = 3 , byrow = T)
severity <- classify (dnbr, rclas)
SCcolors = c ("#008080" ,
"#5f9ea0" ,
"#e0e0e0" ,
"#a0522d" ,
"#8b0000" )
SCnames = c ("Unburned" ,
"Low" ,
"Moderate" ,
"High" ,
"> Green" )
severity_df <- rasterdf (severity)
ggplot (severity_df) +
geom_raster (aes (x = x,
y = y,
fill = as.character (value))) +
scale_fill_manual (name = "Severity Class" ,
values = SCcolors,
labels = SCnames,
na.translate = FALSE ) +
annotation_scale (location = 'bl' ) +
coord_fixed (expand = F) +
theme_void ()
Code
severity_counts <- severity_df %>%
count (value) %>%
mutate (
percentage = (n / sum (n)) * 100 , # Convert counts to percentages
value = as.character (value) # Ensure value is a character for plotting
)
# Create the bar chart showing percentages
ggplot (severity_counts, aes (x = value, y = percentage, fill = value)) +
geom_bar (stat = "identity" ) +
scale_fill_manual (name = "Severity Class" , values = SCcolors, labels = SCnames) +
labs (x = "Severity Class" ,
y = "Percentage of Pixels" ,
title = "dNBR Percentages by Severity Class" ,
subtitle = "A0103: Paso Caballos" ) +
theme_minimal () +
geom_text (aes (label = sprintf ("%.1f%%" , percentage)), # Format label to 1 decimal place
vjust = 1 , size = 5 ) +
theme (
plot.title = element_text (hjust = 0.5 , size = 18 , face = "bold" ), # Center title & increase size
axis.title = element_text (size = 14 )
)
Code
severity_counts <- severity_df %>%
count (value) %>%
mutate (
percentage = (n / sum (n)) * 100 , # Convert counts to percentages
value = as.character (value) # Ensure value is a character for plotting
) |>
filter (! is.na (value))
# Create the bar chart showing percentages
severity_plot <- ggplot (severity_counts, aes (x = value, y = percentage, fill = value)) +
geom_bar (stat = "identity" ) +
scale_fill_manual (name = "Severity Class" , values = SCcolors, labels = SCnames) +
labs (x = "Severity Class" ,
y = "Percentage of Pixels" ,
title = "dNBR Percentages by Severity Class" ,
subtitle = "A0105: Sierra del Lacandon" ) +
theme_minimal () +
geom_text (aes (label = sprintf ("%.1f%%" , percentage)), # Format label to 1 decimal place
vjust = 1 , size = 5 ) +
theme (
plot.title = element_text (hjust = 0.5 , size = 18 , face = "bold" ), # Center title & increase size
plot.subtitle = element_text (hjust = 0.5 , size = 14 )
)
ggsave ("outputs/A0105-dNBR-Percentages-by-Severity-Class.jpg" ,
plot = severity_plot, width = 8 , height = 6 , dpi = 300 )
Code
# # Define severity classification matrix
# rclas <- matrix(c(-Inf, -970, NA, # Missing data
# -970, -100, 5, # Increased greenness
# -100, 80, 1, # Unburned
# 80, 265, 2, # Low severity
# 265, 490, 3, # Moderate severity
# 490, Inf, 4), # High severity
# ncol = 3, byrow = TRUE)
# Color palette
severity_colors = c ("#008080" ,
"#5f9ea0" ,
"#e0e0e0" ,
"#a0522d" ,
"#8b0000" )
severity_names = c ("Unburned" ,
"Low" ,
"Moderate" ,
"High" ,
"> Green" )
# Create leaflet map
leaflet () %>%
addTiles () %>%
addRasterImage (
severity,
colors = severity_colors,
opacity = 0.7 ,
group = "Severity Classes"
) %>%
addLegend (
position = "bottomright" ,
colors = severity_colors,
labels = severity_names,
title = "Severity Class"
)